class'Polynom'
function Polynom:__init(...)
	self.coefficients = table.pack(...)
	self:FixCoefficients()
	return self
end

function Polynom:FixCoefficients()
	while self.coefficients[1]==0 do 
		table.remove(self.coefficients,1) 
	end
end

function Polynom:GetValue(n)
	local nval, res = #self.coefficients, nil
	if n and nval>0 then
		if math.abs(n)==math.huge then --Limes
			res = self.coefficients[1]*n^(nval-1) or nil
		else
			for i=1, nval do
				res = (res or 0) + self.coefficients[i]*n^(nval-i)
			end
		end
	end
	return res
end

function Polynom:Derived()
	local coeff = {}
	local nval = #self.coefficients
	for i=1, nval-1 do
		coeff[i] = (nval-i)*self.coefficients[i]
	end
	return nval>0 and Polynom(table.unpack(coeff)) or nil
end

function Polynom:Integrated(c)
	local coeff = {}
	local nval = #self.coefficients
	for i=1, nval do
		coeff[i] = self.coefficients[i]/(nval-i+1)
	end
	coeff[nval+1] = c or 0
	return Polynom(table.unpack(coeff))
end

function Polynom:GetInflections(derivate2)
	derivate2 = derivate2 or self:Derived():Derived()
	return derivate2:Zeros()
end

function Polynom:GetLocalExtrema(derivate)
	derivate = derivate or self:Derived()
	local extrema = derivate:Zeros()
	table.sort(extrema)
	return extrema
end

function Polynom:GetLocalMinima(derivate, derivate2)
	derivate = derivate or self:Derived() 
	derivate2 = derivate2 or derivate:Derived()
	local res = {} 
	for i, v in ipairs(self:GetLocalExtrema(derivate) or {}) do
		res[#res+1] = derivate2:GetValue(v)>0 and v or nil
	end
	return res
end

function Polynom:GetLocalMaxima(derivate, derivate2)
	derivate = derivate or self:Derived() 
	derivate2 = derivate2 or derivate:Derived()
	local res = {} 
	for i, v in ipairs(self:GetLocalExtrema(derivate) or {}) do
		res[#res+1] = derivate2:GetValue(v)<0 and v or nil
	end
	return res
end

function Polynom:GetGlobalExtrema(derivate)
	local compv = {
		{-math.huge, self:GetValue(-math.huge)},
		{math.huge, self:GetValue(math.huge)}
	}
	local res = {}
	for i, v in pairs(self:GetLocalExtrema(derivate)) do
		compv[#compv+1] = {v,self:GetValue(v)}
	end
	table.sort(compv,function(a,b)
		return a[2]<b[2]
	end)
	local miny, maxy = compv[1], compv[#compv]
	for i=1, #compv do
		local x, y = table.unpack(compv[i])
		if y==miny or y==maxy then
			res[#res+1] = x
		end
	end
	return res 
end

function Polynom:GetGlobalMinimum(derivate, derivate2)
	local x1, x2, x3 = -math.huge, math.min(math.huge,table.unpack(self:GetLocalMinima(derivate, derivate2) or {})), math.huge
	local y1, y2, y3 = self:GetValue(x1), self:GetValue(x2), self:GetValue(x3)
	local xR, yR = x1, y1
	if y2<yR then
		xR, yR = x2, y2
	end
	if y3<yR then
		xR, yR = x3, y3
	end
	return xR, yR
end

function Polynom:GetGlobalMaximum(derivate, derivate2)
	local x1, x2, x3 = -math.huge, math.max(-math.huge,table.unpack(self:GetLocalMaxima(derivate, derivate2) or {})), math.huge
	local y1, y2, y3 = self:GetValue(x1), self:GetValue(x2), self:GetValue(x3)
	local xR, yR = x1, y1
	if y2>yR then
		xR, yR = x2, y2
	end
	if y3>yR then
		xR, yR = x3, y3
	end
	return xR, yR
end

function Polynom:RealRoots(derivate)
	local res, m = {}, #self.coefficients
	self:FixCoefficients()
	if m==2 then
		local m, c = self.coefficients[1], self.coefficients[2]
		res[1] = -c/m
	elseif m==3 then
		local a, b, c = self.coefficients[1], self.coefficients[2], self.coefficients[3]
		local sqr = b*b-4*a*c
		if sqr>=0 then
			local nom, denom = math.sqrt(sqr), 2*a
			res[1] = (nom-b)/denom
			res[2] = (-nom-b)/denom or nil
		end
	elseif m>3 then --Aproximation using Newtons approximation
		derivate = derivate or self:Derived()
		local epsilon = 2^-23 --float mantissa (precision near 0) is 2^24-1
		local function approximate(p1, p2, xMin, xMax)
			local x0 = p1.y~=p2.y and (p2.x*p1.y-p1.x*p2.y)/(p1.y-p2.y) or nil
			if x0 and xMin<=x0 and x0<=xMax then
				local y, oldVal, newVal, i = self:GetValue(x0), nil, x0, 1
				while math.abs(y or math.huge)>epsilon do
					if i>100 then break OutputDebugString("Newton Approximation Limit") end
					oldVal = newVal
					y = self:GetValue(oldVal)
					newVal = oldVal - y/derivate:GetValue(oldVal)
					i = i + 1
				end
				return newVal
			end
		end
		local extrema = self:GetLocalExtrema(derivate)
		local firstExtrema, lastExtrema = extrema[1], extrema[#extrema]
		local theta = firstExtrema~=lastExtrema and math.abs(lastExtrema-firstExtrema) or 1/(math.pi*3)
		local p1, p2 = {x = firstExtrema-theta, y = self:GetValue(firstExtrema-theta)}, {x = firstExtrema, y = self:GetValue(firstExtrema)}
		while (p1.y==p2.y) do
			p1 = {x = firstExtrema-theta, y = self:GetValue(firstExtrema-theta)}
			theta = theta*2
		end
		res[#res+1] = approximate(p1, p2, -math.huge, p2.x)
		p1 = p2
		for i=2, #extrema do
			p2 = {x = extrema[i], y = self:GetValue(extrema[i])}
			res[#res+1] = approximate(p1, p2, p1.x, p2.x)
			p1 = p2
		end
		while (p1.y==p2.y) do
			p2 = {x = lastExtrema+theta, y = self:GetValue(lastExtrema+theta)}
			theta = theta*theta
		end
		res[#res+1] = approximate(p1, p2, p1.x, math.huge)
	end
	return res
end

function Polynom:Zeros(derivate)
	local t, roots, lastRoot = {}, self:RealRoots(derivate), nil
	for i=1, #roots do
		local v = roots[i]
		if not lastRoot or math.abs(lastRoot-v)>2^-23 then
			t[#t+1] = v
			lastRoot = v
		end
	end
	return t
end